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1. Abstract 

In this paper the Navier-Stokes-a (NS-a) model is considered within a large-eddy 
simulation framework. An investigation is carried out using fully-developed turbu- 
lent channel flow at a fairly low Reynolds number. This is a flow where diffusion 
plays a prominent role, and presents a challenge to the nonlinear model investi- 
gated here. It is found that when a\ is based on the mesh spacing, the NS-a model 
has a tendency to tilt spanwise vorticity in the streamwise direction, leading to 
high skin friction. This is due to interaction between the spanwise vorticity, the 
model, and the streamwise streaks. To overcome this problem a\ is damped in the 
streak affected region. Results overall demonstrate the potential of the model to 
reproduce some features of the DNS (helicity statistics and small-scale features), 
but more work is required before the full potential of the model can be achieved. 
In addition to the channel flow investigation, a derivation of the governing equa- 
tions using Hamilton's principle is given. The derivation is intended to be clear 
and accessible to a wide audience, and contains a new interpretation of the model 
parameter. 

Keywords: Navier-Stokes-alpha; large-eddy simulation; plane channel flow; regularization 
model; subgrid-scale model 



2. Introduction 

Traditionally, turbulence models are derived by applying averaging (RANS) 
or filtering (LES) techniques to the Navier-Stokes equations. This results in a 
momentum equation with an unclosed term, known as the Reynolds stress or 
subgrid stress, that must be modeled. Numerous models have been proposed over 
the years il-Si]. The majority of these models employ an eddy viscosity ansatz. 
This is well-founded in the sense that an eddy viscosity is a reasonable model for 
the energy drain provided by the small scales that have been removed during the 
filtering or averaging procedure, and is popular in part because adding viscosity 
generally renders a simulation more stable. However, the shortcomings of the eddy 
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viscosity approach are well known. Linear eddy viscosity models employ a simple 
constitutive relationship that assumes alignment between the subgrid stress, Tjj, 
and the strain rate. It has been found both in analysis of DNS data and in 
experimental studies that this is far from the truth [4]. Eddy viscosity models are 
also strictly dissipative for positive viscosities, and unstable for negative ones. This 
means they cannot capture the reverse energy transfer from small to large scales, 
known as backscatter. Although energy transfer is on average from large to small 
scales in three-dimensional turbulence, there are a number of flows where local 
backscatter effects are important. Examples from shear flows include the later 
stages of boundary layer transition 0], hairpin vortices in the near- wall region 

0, Q] and vortex pairing in mixing layers [1]. A popular method to incorporate 
backscatter is by adding a stochastic forcing term to the eddy viscosity model 

01. The rationale behind this is that while a dissipative model can capture the 
mean forward transfer, the subgrid stress exhibits significant fluctuations about 
this mean, and it is these fluctuations that are responsible for the backscatter. 
In practice though, backscatter tends to be strongly correlated with coherent 
structures, leading some to hypothesize that a deterministic model may be more 
appropriate 



One flow where traditional eddy viscosity models have difficulty is turbulent 
channel flow. In this case, the eddy viscosity needs to be reduced close to the wall, 
or a dynamic procedure needs to be used, to avoid damping out the turbulence. In 
this paper we study the turbulent channel flow using the NS-a model. The NS-a 
model is different from an eddy viscosity model in that instead of adding an eddy 
viscosity term to a filtered momentum equation, it is a nonlinear regularization. 
The model can be thought of most intuitively as a vorticity regularization. The 
vorticity equation for the inviscid form of the NS-a model is (in this paper the 
use of repeated indices implies a summation, unless otherwise stated), 

dt dxj dxj ' 

In this equation the background fiow is smoothed, which means the velocity 
gradients become less effective at stretching and tilting the vortices. In turn this 
suppresses the forward transfer of energy and prevents the creation of smaller and 
smaller scales, hence eliminating the need to model the effects of these scales when 
we carry out a coarse grid numerical simulation. In this way additional viscosity 
is not needed, per se. This phenomenological view of the model is supported by 
Fourier transform analysis of the nonlinear terms l3, U | , which shows that if the 



underlying dynamics follow those of the Navier-Stokes equations, then the system 
stops transferring energy to small scales when a certain wavenumber, say /cq, is 
exceeded. At the same time, it does not stop the backscatter. Thus we do not 
need to model the missing backscatter because it has not been removed. 

In spite of its intuitive application ctS cl turbulence model [i3-[ill , there 
have been only a few attempts in the literature to use the NS-a model outside of 
idealized box turbulence experiments [H, [iMl. Geurts and Holm ^ used the 



model to capture temporal transition in a mixing layer, while Holm and Nadiga 
(21I were able to produce a four-gyre structure on a coarse mesh that would only 
produce two gyres when viscosity was used as the closure. Recently, the NS-a 
model was incorporated into a primitive-equation ocean model where it was found 
to produce energetic eddies at a coarse resolution where eddy viscosity approaches 
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failed ^22] , In all of these studies they maintained a constant value for the model 
parameter a, as a fraction of the mesh spacing. This was done because a relates 
the smoothed and unsmoothed velocities and can be interpreted as a filter width. 
However, we expect in flows which are highly anisotropic, such as near wall flows, 
that we will not be able to maintain a constant value of a. This topic was explored 



by Scott and Lien [231] the NS-a model was applied to a lid-driven cavity flow, 
and both a as a function of the mesh spacing and a flow dependent version were 
tested. The flow dependent version was found to reproduce the results of the DNS 
fairly well. 

Another study that used a non-constant a was an investigation of turbulent 
channel flow by Zhao and Mohseni [2^, where a dynamic version of the NS-ct 
model was developed and tested in an a priori manner, in which was cal- 
culated, but was not fed back to the flow. A later investigation of the channel 
flow where the model was tested a posteriori found the model to produce high 
spanwise fluctuations [25]. One of the objectives of the present study was to 
determine the source of this problem. We will show here in section 14.4.21 that 
this occurs because the model has a bias towards tilting vorticity in the stream- 
wise direction close to solid walls. Our response to this bias is given in section 14.4.31 

Before this, a derivation of the governing equations is presented in Section 
[3l followed by the formulation of the subgrid model and description of the test 
case. The derivation is comprised in part of principles seen before [121 . 

M,^, and 



is intended to be more accessible to the non-mathematician than that commonly 
encountered in the literature Unlike previously published derivations. 



the current one shows explicitly the steps used in varying the action, and we hope 
it will serve as a useful basis for further extensions of the model to different flows. 
The derivation also contains a new interpretation of the model parameter, showing 
how it can be related to a particle displacement error, and shows explicitly how 
the variation with respect to this parameter is carried out. 



3. Derivation 

The NS-a equations differ from other approaches to turbulence modelling in that 
the effects of turbulence are introduced at the level of the variational principle. This 
is done here using Hamilton's principle, which is a variational principle that leads to 
Newton's second law. To incorporate the effects of turbulence within the framework 
of Hamilton's principle, consider that in the material description of a fluid, the state 
of a fluid particle with label a is specified by the particle displacement x(a, t), and 
velocity, x(a, t). A momentum equation arises when the first variation of the action 

S[x,x\ = [ [ lm{x,x)(fadt, (2) 

Jti Jv{a) 

is set to zero. In ([2|) is the material representation of the Lagrangian density 
(Lagrangian/unit volume), and is the difference between the kinetic energy T and 
the potential energy <I>. Turbulence can be incorporated within this framework by 
adding a random component to the displacement of a fluid particle, and, given 
that this random component is a function of time, to its velocity. This is the 
method pursued by Marsden and Shkoller [2§] and Holm . 
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In the following we will pursue a different approach and work with the Eu- 
lerian description of a fluid, where the state of the (isentropic) fluid is described 
by u(x,t) and p(x,t). This is the description that is usually used in develop- 
ing a turbulence model. For example, the well-known RANS equations for an 
incompressible flow are developed by decomposing the velocity field into mean 
and fluctuating components, substituting into the Navier-Stokes equations, and 
averaging. To facilitate the derivation of the model we have split it into the five 
sections. Here a brief description of each section is given. 

• Definition of the Lagrangian In this section we describe the particular form of 
the Lagrangian we are using, which is for constant density, incompressible flow 
with no sources of potential energy. To apply the model to different flows, the 
definition of the Lagrangian must be modified. Here we write the Lagrangian in 
Eulerian coordinates, but it is also possible to work with material coordinates. 

• Incorporation of turbulence into the Lagrangian In this section we describe the 
definition of the velocity fluctuation. This is the only approximation that enters 
the NS-a model. This section follows that given in Holm [13], with a different 
interpretation of the model parameter. 

• Varying the action In this section the first variation of the action is taken. This 
is an application of the calculus of variations. If different boundary conditions 
are used for the model parameter, this part should be modified. 

• Definition of the variations of the Eulerian coordinates In this section the rela- 
tionship between a trajectory variation and the Eulerian variables are used to 
define their variations. This section is particular to the use of Eulerian coordi- 
nates. 

• Setting the variation to zero. When the first variation is set to zero, we arrive at 
the momentum equation. 

3.1. Definition of the Lagrangian 

In Eulerian coordinates the action principle is (c.f. [H, [sij), 

S[u,p,s]= / / l{u, p, s)d^x dt. (3) 
Jti Jv(x) 

with Lagrangian density 

/ = ^^^Ui{x, t)ui{x, t) - E{p{x, t), s{x, t)) - V{x, t) (4) 

where E is the internal energy, s is the entropy and (j) is the potential energy. Here 
we will consider an incompressible, constant density fluid of uniform entropy with 
no sources of potential energy. In the action principle we then remove the internal 
and potential energy functions and add an equation constraining the density to be 
constant 

^P = l J^(^^^Y^Ui{x,t)u,{x,t)+p{po-p{x,t))^d^xdt. (5) 

Here p is the pressure, and is a Lagrangian multiplier. Using conservation of mass 
we can relate the density ratio to the volume element D by D = p/po, where po is 
the reference density. The volume element is defined as the ratio of the volume 
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in the initial configuration to that in the current configuration 

D = det{da/dx) (6) 
Using this we arrive at the action 

where we have divided through by po- 

3.2. Incorporation of turbulence into the Lagrangian 

To incorporate turbulence the velocity is then expressed as the sum of a mean 
component and a random fluctuation, in a similar manner as what is done in 
RANS (here (/> is a random variable, or a fast time-scale [l^ ) 

Ui{x,t](l)) = Ui{x,t) + u'i{x,t]4)), (8) 

where the averaging operator (•) and u are defined as [l^ 

1 [T 

u{x,t) = {u{x,t;uj)) = lim — / u{x , t; u) duj . (9) 

T^oo T Jq 

The only approximation in the NS-a model comes in the definition of the velocity 
fluctuation. For example, we can write (to first order) 

iii{x + ^) = Ui{x) + Cj-^. (10) 
dxj 

Defining the velocity fluctuation as the difference between our averaged velocity at 
two points, u{x) and u{x + ^), gives 

, . dui 

For example, we can consider that if we are sitting at a field point x occupied 
by a particle with velocity u{x) at time t, and then at a later time our field 
point is occupied by a particle that was previously at x + ^, and has velocity 
u{x + ^), the velocity fiuctuation is then given by (jlip . The same expression for 
the velocity fluctuation can be derived in a similar manner by expanding Eulerian 
and Lagrangian velocities in terms of ^ [13] • 

An alternative interpretation of this picture can be found by looking at ^ 
as the error between true and modelled trajectories. Note that in this discussion 
on the interpretation of the model parameter as an error there is no summation 
on repeated indices. Given that 

Dx , , 
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a simple first order discretization of the true trajectory Xt and the modehed tra- 
jectory Xm would be 



t \ t I . (13) 

™n n-1 I n-1 fa-n-l^ /\ + 

where the superscript n indicates the time level. Defining the error as the diff'erence 



between the true and modehed trajectory we find [32l | 



e = + {ur'i^r') - <~H<-')) a*. (u) 

We can relate the true velocity to that at the modelled particle location using the 
definition of the error 

<-i«-i) = <~i«-i+r-') 

(15) 

«<-i(x«-i)+r-'-v<-^ 

Then, split the true velocity into a large and small scale component (i.e. Ut = 
ui + Us), and assume the large component is equal to the modeled field to obtain 
(neglecting products of $, and the small scale velocity) 

^ ^n-1 ^ ^^n-l(^n-l) ^ ^n-1 . Vw^-^) At. (16) 

This is a discrete form of the following equation, where we assume that ui u 
and Us — u' 

^ = «' + ^V^Z. (17) 

Setting D^/Dt = 0, which means the error is frozen along a particle trajectory, 
or that all of the fluctuation is contained in the Eulerian field we arrive at 
the deflnition of the velocity fluctuation from before, equation (llip . Decomposing 
the velocity in the Lagrangian in equation ([7]) into a mean and a fluctuation, 
substituting the velocity fluctuation from (lllh . and averaging (using {$,) = 0, which 
means the error is unbiased) yields the averaged Lagrangian, 

A. (18) 



Here we have followed the notation used in Holm where the averaged velocity is 
u, but where we keep the brackets for the averaged displacement covariance {^kCi)- 
We have in equation (jlSh essentially our turbulence model. We can see that the 
kinetic energy is composed of two parts, the first is the kinetic energy of the mean 
fiow, and the second part is the kinetic energy of the fiuctuating component (or 
eddy kinetic energy). If we consider only the diagonal components of {(,k(,i), the 
energy due to the fiuctuating part will remain positive. 
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3.3. Varying the action 

To obtain our momentum equation we need to set the first variation of the action 
to zero. The action is defined in the usual manner 

S= [\L)dt, (19) 
Jti 

with (L) defined in equation (|18p . The first variation of the action is 

where (/) is the Lagrangian density (Lagrangian/unit volume). The partial deriva- 
tives with respect to the volume element and particle displacement are 

d{l) ^ UjUj {^kCi) duj duj 

dD ~ 2 2 dxkdxi ^' ^ ^ 

d{l) D dui du. 



d{(.k^i) 2 dxk dxi ' 
For the velocity 



(22) 



/ -^OUi dV = / Duid Ui + D - — [dUi) h t; — t;— [.OUi) d-'x 

Jv 9ui Jv 2 \dxk dxi dxkoxi J 



(23) 



: DuMi + D{ikii) (1^^ (-^^O) (24) 



(25) 

where we have applied symmetry of the particle displacement covariance {£,k(,i), 
and applied integration by parts in the last step. If we apply either a periodic 
boundary condition or the constraint that the normal component of {^k^i) is zero, 
the surface integral in (j25p is zero and we have 



3.4. Definition of the variations of Eulerian coordinates 

We have now defined the partial derivatives which appear in our varied action 
(I20p , it still remains to define the variations of the Eulerian coordinates 6ui , 6D 
and 5{^k^i). If we were working with the material representation, given by equation 
([2]), variations would have been taken with respect to the particle position, and 
all that would remain to be done would be to set the first variation to zero. In 
the spatial picture we have followed here variations of the Eulerian coordinates 
{u, p, s) are taken at a fixed point, although the final goal is the same as in the 
material representation, to find the trajectory for which the action is stationary. 
The first question is then, how are variations of a particle trajectory reflected in 
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the Eulerian coordinates? 

To connect the two we define a function that relates the field position x to 
a label a. Such a function is given by the trajectory, which we will denote by rj. 
More formally, let rj be the function that maps particles with labels a to the field 
points they occupy at time t. For our purposes here we will assume this map is 
one-to-one, invertible and sufficiently smooth that we may differentiate it as many 
times as necessary. The particle position x is 

x = ri{a,t). (27) 

Similarly, let r]~^ be the map that tells you the label a of the particle occupying 
the field point x at time t, 

a = r]-^{x,t). (28) 

The Eulerian and Lagrangian velocities at a given point are related through the 
identity 

u = ^{r]-\x,t),t), (29) 



u = i] or] ^ (30) 



or 



where the o operator denotes a composition of maps [33], and the dot indicates a 
time derivative. Consider that if we only know the mapping rj, then to find the 
velocity, u{x, t), at a given field point x, we can evaluate a = r]^^{x, t) at our field 
point to find the particle occupying that point at time t. Knowing the particle 
(denoted by the label a) we can then evaluate rj{a, t) to get the velocity at that 
point. The central idea to keep in mind here is that in the Eulerian framework 
the dependence of the velocity field on the particle trajectory comes into play in 
two places, one in calculating the rate of change of the trajectory and the other in 
evaluating the label. This means when we vary the velocity field we need to take 
both of these into account. 

The trajectory is related to the volume element according to 



Finally, the displacement covariance equation. 



Dt 

can be written in material form as 



0, (32) 



mi) ° V = mi)o (33) 



where the subscript o denotes the initial value. This equation tells us the 
displacement covariance at the current particle location r/ is equal to the initial 
displacement covariance, in other words that the displacement covariance is 
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preserved along trajectories. From these expressions for the velocity (I30p . volume 
element (f3T]l , and displacement covariance (f33l) we will be able to relate variations 
of the Eulerian coordinates to the particle trajectory. 

As an example consider the variation of the displacement covariance. We 
start with the definition of the particle trajectory variation |27,] 



6r] :- 



d 
Te 



{T] + e5ri). 



(34) 



e=0 



Variations of the Eulerian coordinates, u,D and {S,k(,i), are defined in a similar 
way. Now vary both sides of (jSSp and differentiate with respect to e noting that 
the RHS is constant 



d 
Te 



e=0 



Denoting the varied quantities as 

r]'^ = T] + eSrj 
and carrying out the differentiation gives 

V(6ez) • + 5{ikii) o 77 = 0. 
Composing both sides with r/~^ 

5{ikii) = -y{ikii) ■ w 

where w is the trajectory variation expressed at a field point, 

w = Srj o T]~^ . 



(35) 



(36) 
(37) 



(38) 



(39) 



(40) 



In a similar manner it can be shown that the velocity and volume element variations 
are 0,[27| 



„ dw _ ^ ^_ 
ou = — — h u ■ \/w — w ■ \/u, 
ot 



(41) 



5D = -V ■ {Dw) . (42) 

The first two terms in the velocity variation are due to the rate of change of the 
trajectory variation, with the label fixed, while the last term is due to the variation 
with respect to the label . This is what we meant earlier when we said earlier 

that both the trajectory and label need to be varied. Expanding the variation of the 
volume element gives two terms that can be interpreted similarly. A more detailed 
discussion on the variation of Eulerian quantities can be found in Bretherton [26| . 
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3.5. Setting the first variation to zero 

To proceed with setting the first variation of the action, given by equation ()20p to 
zero, we substitute the variations given by equations (j4ip . (j42p . (j39p into the varied 
action (f20|l . 



6S 



d(l) ( dwi _ dwi 

' + Uj 



duj \ dt 



Integrating by parts and changing sign 



d_ 

dxi 



d{l) d{l) d{^k(i) 



dD di^kW dxi 



(43) 



Wj I d^x dt. 



6S 



h JV 

Jv dui 



d_ _ _d_\ a(0 d{l) duj 
dt dxj J dui duj dxi 



■Wi dV 



-^—UjWi dAj + 



dxi dD 5(66) dxi 

*^ d{i) 



-Dwj dAj 



II 



III 



Wjd^x dt 



(44) 



The last three terms are zero for the following reasons: 

(I) variations are zero at beginning and end times (same as for Newton's law), 

(II) velocity is either periodic or has zero normal component for a solid surface, 

(III) trajectory variation {w) is tangent to the bounding surface 34l |. 



Substituting in the partial derivatives ([^ . (piD . ([2^ into taking the 

limit as d^x,Wi, and dt go to zero (c.f. Gelfand [35.]) and applying the constraint 
D = 1, yields the momentum equation 



dt 



— - 

dxi dxj 



dvi _ dv, 



dp" _ 1 5(66) dum dum 
dxi 2 dxi dxk dxi 



(45) 



where the following variables have been defined 



Vi = Ui- —— {Ck(.l)^ 

dxk V 



p"" =p - ^UmUm - (66 



dui dui 
dxu dxi 



(46) 
(47) 



Note that there are two velocities, related through the Helmholtz operator. 



Vi 



— f (66)— 

dxk \ ' dxi 



Ui 



(48) 



H 



If we impose isotropy (66) = ct^^ki s^-nd assume is constant, we arrive at the 
inviscid form of the NS-a equations found in the literature 20|, ]21 \ 



dvj 
dt 



+ u. 



dxi 



+ Vj 



9uj 
dxi 



dxi 



(49) 



June 13, 2010 17;26 Journal of Turbulence hpchan5 



12 

with 



dxl 



In Fourier space 



k{k) = T^^, (51) 
1 + fer 



from which it is clear that the smoothed velocity is low-pass filtered since the 
high- wave number components are attenuated. Prom this relationship we can see 
that a can be interpreted as a filter width. 

The continuity equation does not come from the variational principle, but 
from taking the time derivative of the volume element [l^ 



dD 

— + V • {Du) = 0. (52) 

Imposing the constraint D = 1 then gives V-w = 0, orV-'U = when we 
recognize u is the smoothed velocity according to the Helmholtz operation. 

There are a few things to note here. The first is that the averaged velocity 
u becomes the smoothed velocity when we consider that v and u are related 
through a Helmholtz operator. This is what is meant in the literature by 'temporal 
averaging in the variational principle implies a spatial smoothing in the momentum 
equation' [l^, although one could anticipate this from the expression for the 
velocity fiuctuation, which was derived using a spatial Taylor series expansion (jlOp . 
Another aspect to note is that if we had not considered the functional dependence 
on {^kCi) we would not have obtained the d{^k^i)/dxi term in the final momentum 
equation (j45p . This term is necessary to conserve momentum, which you can see 
either by considering Noether's theorem for the action principle, or by removing 
this term and trying to write the momentum equation in conservative form. 
Some studies of the NS-a equations have not included this term in their analyses 



yielding incorrect results, as pointed out in the literature [36|, |33]. Finally, by 
following through with this method we are able to understand how the boundary 
conditions for the NS-a equations arise. 

Note that other methods can also be used to derive the equations from 
Hamilton's principle [i3, m, 0, [si]. The method given here is both straightfor- 



ward and general enough to allow the model parameter {^k^i) to be non-constant 
and anisotropic. Extensions of the model to stratified and rotating fiows are 
given by Holm [H, [3§], and can be obtained using the methods used here. For 
alternative examples on the use of variational principles in fluid mechanics see 
Salmon [sH and Finlayson [40]. When using their methods the advection equation, 
D{^k^i) / Dt = 0, needs to be added as a constraint equation to obtain the last 
term on the RHS of Equation (I45p that contains the gradient of the particle 
displacement covariance. 
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4. Demonstration 

4.1. Model Formulation 

The formulation of the NS-a model has been described in [13] and is briefly re- 
viewed here for continuity. To investigate the NS-a model numerically we work 
within an LES template and develop an equation with the smoothed velocity as 
the dependent variable. To do this, first add a viscous term to (j45p and then rewrite 
the equation in momentum-conservation form [iS?] (now replacing u with u) 

dvj ^ ^ dvj _ dp ^ d ^ dumdum \ ^ ^ p'^Vj 
dt dxj dxi dxj \ dxi dxk J dx\ ' 

Then, to write the substantial derivative entirely in terms of the smoothed velocity, 
rewrite the advective terms in (1531) as, 



Here [D/Dt, H] is the commutator between the material derivative and the 
Helmholtz operator, H from equation (j48p . 

[D/Dt, H]ui = D/Dt{H{ui)) - H{D/Dt{ui)), (55) 

where H{ui) = Vi. Note that the substantial derivative is defined with the smoothed 
velocity, D/Dt = dt + Ujdj. The momentum equation can then be written as, 

^j^u — = H-'^ (-^ + — ( {^,^^.)^^^^\ - [D/Dt H]u^ 

dt '' dxj \ dxi dxj \ •' dxi dxk J dxi ' / ' 

(56) 

Expanding the commutator, applying D{^k£_i)/Dt = 0, the momentum equation 
can then be written 

dui ^ _ dui dp ^ d'^Ui ( dm. 



_l_ _|_ I' _ j j (57) 

dt dxi dxi dxi \ dxi ' 



The subgrid stress is 



dui duj / f > \ duk dui , ^ ^ > dum dum 



Instead of using the full anisotropic model, in the following demonstration we con- 
sider a simplified version that arises when only the diagonal components of {£,kS,i) 
are retained, which ensures that the kinetic energy in the Lagrangian (equation 
(fT8]) ) is positive. Denoting a| = (CkCk) we arrive at our subgrid stress 

2^ dUidUj 2e duk dui dumdurn 



^ij 

This is the subgrid stress that would result if the Helmholtz operator is considered 
as being equivalent to the composition of three one-dimensional, symmetric filters 
(the off-diagonal components now being zero). A similar filter has been used in 
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the Tensor-Diffusivity model [4l|- Here it reduces the cost of the model such that, 
when the explicit filter is applied by solving the Helmholtz equation using Fourier 
transforms, the model adds approximately 30% to the total computational time, 
similar to what is found in other studies that used a constant, isotropic model 
parameter, a'^Ski = (66) 



4.2. Description of the test Case 

The test case chosen here is turbulent channel flow. The focus here is on wall- 
resolving LES, thus we are going to consider Reynolds numbers at the low end of 
the turbulent regime. This is a challenging test case for the present model because 
it is a model with a modified nonlinearity, while in this test case diffusion plays a 
prominent role. Channel fiow with Reynolds number Rcr = 180 is studied here 
using a second-order finite volume method (43| . Periodic boundary conditions 
are applied in the homogeneous directions (streamwise and spanwise) while 
no-slip conditions are used for the solid boundaries located at y = ±H, where 
H is the channel half- height. The mesh is uniformly spaced in the homogeneous 
directions but stretched in the wall-normal directions using a hyperbolic tangent 
profile. To enable a variety of subgrid and numerical resolutions to be tested, 
the investigation was done primarily using the minimal channel flow [3]. This is 
the smallest domain for which the near wall cycle is able to sustain turbulence. 
The near wall cycle consists of interactions between low-speed spanwise streaks. 



streamwise vortices, and hairpin vortices |45l |. To sustain turbulence, the channel 
must be wide enough to contain a low-speed streak. For this purpose we chose a 
channel of dimension (vr, 2, O.Stt), where the non-dimensionalization is with respect 
to the channel half height. Various mesh resolutions were tested, summarized in 
Table [TJ In our simulations a constant mean mass flux was enforced at a Reynolds 
number of 4160 based on the centerline velocity of a laminar flow and the channel 
half-height H. This is equivalent to a bulk flow Reynolds number of Rcb = 2773 
or Rcr = 180. 

The flow was initialized using a parabolic profile with a superimposed Tollmien- 
Schlicting (T-S) wave to provide a 2D disturbance. A T-S wave with amplitude 
of 10% of the centerline velocity and wavenumber 2 (made dimensionless with 
the channel half-height) was found to bring the flow to a turbulent state quickly. 
This method was preferred over white noise because it was found the noise had 
a tendency to require a longer time to reach a turbulent state. By initializing 
the flow with a large scale disturbance, nonlinear interactions quickly generate a 
cascade of energy towards the small scales. 

In the results an averaged quantity is denoted by an overbar, and a fluctu- 
ation about this state is denoted with a prime. For the velocity, vorticity and other 
profiles reported (quantities that are a function of the vertical coordinate) the 
averaging is taken over the statistically homogeneous directions x and z as well as 
with time to increase the statistical sample. Quantities are non-dimensionalized 
using the channel half-height, H ^ and the shear velocity u^- = \l^w , where is the 



wall shear stress = v 



These non-dimensional quantities are n"*" = u/ur, 
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4.3. Definition of aj^ 

To specify a\ here 
based a\ on the mesh spacing 



15 

2 



To specify a\ here as a first step we followed a conventional LES approach and 



C {hi) (60) 



where hk is the grid spacing in the fc-direction and C is a constant denoting what 
fraction of the grid spacing to use. Because a| can be related to the width, A^, 
of a box filter via a| = A|/24 we choose C = 1/6, which corresponds to 
a filter width which is twice the grid size. To avoid the commutation error that 
arises when a filter with non-uniform widths is used, we chose not to filter in the 
wall-normal direction, and for this reason Oy was set to zero. The filter was applied 
both by solving the Helmholtz equation using Fourier transforms and also by using 



a box filter. Results using the two methods were very similar those shown here 
solve the Helmholtz equation. It should also be noted that we also tried using the 
isotropic model with based on the grid volume, but were not able to achieve 
numerically stable results with that definition of a^. 



4.4. Results 

4-4- 1- Mean flow and energy transfer 

The first quantity that is of interest in turbulent channel flow is the mean flow 
profile, which is related to skin friction. The mean flow profile using the definition 
of a| given in section 14.31 is shown in Figure [T] for the four different meshes listed 
in Table [TJ We can see the velocity is significantly underestimated on all three 
meshes (16, 64, 16; 24, 64, 24; 32, 128, 32) using this definition of a| with C = 1/6. 
This implies that the skin friction is significantly overpredicted. To check the 
robustness of this result we also looked at the effect of refining the mesh while 
keeping the physical size of a| constant. By using a (32, 64, 32) mesh with C = 2/3 
in comparison with a (16,64,16) mesh with C = 1/6, we are able to check the 
effect of increasing the subgrid resolution (the ratio between the filter width and 
the mesh spacing). Geurts and Holm j20l ] found for a temporally evolving mixing 
layer when increasing the subgrid resolution from C = 1/6 to C = 1 they were 
able to reduce the turbulent kinetic energy and bring their simulation results 
into good agreement with DNS data. It can be seen in Figure [T] that in our case 
increasing the subgrid resolution does not improve the mean flow profile. We also 
found the spanwise velocity to be significantly overpredicted close to the wall 
(profiles not shown). A similar result was found in Zhao and Mohseni [? ] in their 
study of a channel fiow. We can see in Figure [2] that the spanwise velocity for the 
NS-a contains more small-scale activity as compared to the case with no model, 
and indicates coherent structures close to the grid scale. This is consistent with 
what is found in other studies [20, 46, 47|. 

Because the NS-a model is non-dissipative, and also contains the backscat- 
ter dynamics, it is possible that if a| is chosen to be too close to the energy 
containing scales this could lead to a build-up of energy. Here we look at the 
subgrid energy transfer term, Tsgs = ^ijdjUi, which represents the energy 
transfer from the resolved to subgrid scales [7|]. Note that rhij = {fnij)-, where 
fhij is defined by equation (j59p . Plots of T^Qg as a function of the wall normal 
distance are shown in Figure [3l Both the total transfer due to the rhij term and 
the individual contributions from the Aij, By and Cij terms (see equation ([59]) ) 
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are shown. Note that the contributions from both the Aij and Bij terms are net 
dissipative, while the Cij term produces net backscatter. For all three terms the 
instantaneous values (not shown) fluctuated about these mean values by an order 
of magnitude, exhibiting both forward transfer and backscatter. 

In Figure [3] for the case where C = 1/6 (filter width of twice the grid size) 
the minimum subgrid transfer (or maximum SGS dissipation) is T^^g ~ —0.06 
which is in good agreement with that reported in the literature from filtering 
DNS data for a channel fiow at the same Reynolds number with the same filter 
@, |4§]. Thus, instead of excessive backscatter or insufficient dissipation, the main 
problem instead is that the dominant physics is too close to the wall. The peak 
transfer in our simulations occurs at y"*" ~ 5, as compared to that in the literature 
at y+ ^ 10. 

4-4- 2- Model bias towards tilting voriticity in the near wall region 

There is a strong correlation between the strength of the streamwise vortices and 
skin friction [49,] . We show here that the vorticity field produced by the NS-a 
model is erroneous, and that this is the cause of the high skin friction (which 
manifests itself as an underpredicted mean fiow profile, shown in Figured]). 

Streamwise and spanwise vorticity fiuctuations, uj'^ = Urmst^/u^^ are shown 
in Figure H for the (24,64,24) and (32,128,32) meshes, both use C = 1/6. The 
minimal channel DNS is in good agreement with the data from the full channel 
DNS [50,] , while the NS-a model significantly overpredicts the streamwise and 
spanwise vorticity fiuctuations very close to the wall. The peak in the streamwise 
vorticity fluctuation at the edge of the buffer layer (at y+ ^ 20 for the DNS) 
is much closer to the wall for the NS-a model (here at ~ 9) and higher in 
magnitude. This peak is indicative of the streamwise vortices in the buffer layer 
[sol ]. According to the streamwise vortex model of Kim et al. j50| , the ratio 
between the streamwise vorticity peak at the wall to that in the buffer layer should 
be ~ 1.3. Here we have instead a ratio of ~ 2 (taking the peak at y~^ ~ 10 to be 
the buffer layer vortices). The wall value of the streamwise vorticity for the NS-a 
model is oj'^lw ~ 0.4, twice that in the DNS where uj^\w ~ 0.18. 

In the Introduction it was highlighted that in the NS-a model vortices are 
tilted and stretched by a smoothed velocity, (see equation ([I])). In the near wall 
region streamwise vorticity is primarily created by tilting of the spanwise vorticity 
into the streamwise direction, through the term in the spanwise vorticity 

equation. Writing the smoothed velocity as Ui = Ui + a\V'^Ui we can see that the 
term responsible for tilting spanwise vorticity into the streamwise direction, w^^, 
is augmented by uizdz [aldzu)- We can compare the NS-a tilting term with that 
from the Navier-Stokes equation by writing the tilting term for the Navier-Stokes 
equation as 



-iy.f. - " (^) , (61) 



NS 



and that for the NS-a equation as 



(62) 
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The most significant contribution to the spanwise velocity gradient, du/dz, 
close to the wall is from the streaks, thus we take this velocity gradient to be 
proportional to the rms streamwise velocity fluctuation divided by the streak 
spacing. The streak spacing can be measured from the two point correlation R 



■uuz 



50(], shown here in Figure [H We found the streak spacing normalized by {v/ur) to 
be narrower with the NS-a model, at ~ 40, than the DNS result by approximately 
a factor of two 0. 

Using the streamwise velocity values from Figure [6l for the minimal chan- 
nel DNS du/dz is approximately 2.7/80 = 0.035, while for the NS-a model it is 
2/40 = 0.05. If we then take the contribution for the second term, dz [a'^d'^^uj, to 
be proportional to {a1/h1) (du/dz) we arrive at the following relationship between 
the two source terms 



^uu. _ ^zidu/dz){l + ayhl) 
<i u^^'idu/dz)^^ ■ 



(63) 



Substituting [dii / dz) / {du / dz)^ ^ ~ 0.05/0.035, al/hl ^ 1/6 and </a;f'^ = 
0.48/0.38 (wall values from Figure U]) into Equation (j63p we arrive at 

= 2.1 (64) 



which agrees well with the values from the minimal channel of 0.40/0.18. 

To investigate the streamwise vortices in the NS-a model further, in Figure 
[7] we compare probability density functions (PDFs) of the streamwise vortex 
inclination angle 6 = a,Tctan{Cjy/Cjx) at two different heights from the wall. For 
the DNS the PDFs were measured using the instantaneous vorticity vector on 
two X — z planes, at vertical locations of ~ 7, and ~ 18. The first is in the 
viscous sublayer, while the second is in the buffer layer. In the viscous sublayer we 
can see two peaks in the PDF at ±90° corresponding to the low and high speed 
streaks. As you move into the buffer region a shoulder appears near 25"/ — 155" 
that corresponds to the streamwise vortices. These results are in good agreement 



with those from the literature |5ll |. 



For the NS-a model to see evidence of streamwise vortices, the PDFs needed to 
be measured closer to the wall, and are shown in Figure [71 At y"*" ~ 5.4 we can see 
there are shoulders near 9 ~ ±90° and more distinct peaks at ~ +20°/ — 160°, 
the former indicating streamwise streaks and the latter indicating streamwise 
vortices. PDFs measured closer to the wall (not shown) had a single peak at ~ 
(zero vertical vorticity, as required by the no-slip condition at the wall). Thus we 
did not see any PDFs indicating a region dominated by low-speed streaks (similar 
to the one at y"*" ~ 7 for the DNS), but instead the region close to the wall shows 
a dominant signature of streamwise vortices. 



^It is not clear at this point why this is the case. One possibility is, given that the streak s pac ing is believed 
to emerge from a secondary instability of the ToUmein-Schlicting wave (Jimenez pg. 219 l4?|), the spacing 
we see here may be related to possible differences that would arise through a stability analysis of the NS-« 
equation as compared to the same anaylysis for the Navier-Stokes equation. For example, it has been shown 
that the model lowers the critical wavenumber for baroclinic instability in a two-layer quasi-geostrophic 
model. Although the initialization here was not representative of a true transition process, there were 
significant differences observed in how the flow became turbulent from the perturbed laminar state when 
the NS-o model was used, as compared to without. 
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4-4-3- Using damping to overcome the model bias 

In the previous section we found that the NS-a model provides an additional 
mechanism for producing atreamwise vorticity in the near-wall region by tilting 
spanwise vorticity directly into the streamwise direction. This is physically 
incorrect in two respects. First, vortex tilting and stretching processes should 
occur farther away from the wall, in the buffer region, not in the viscous sublayer 
which is what we see here. Second, the path of streamwise vorticity creation is 
incorrect. In the literature the streamwise vorticity comes from first lifting up 
the transverse vorticity into the buffer layer and then tilting of the vorticity into 
the streamwise direction [44i], while here we have a direct tilting of spanwise 
vorticity into the streamwise direction. This indicates that either damping 
of in the near wall region or an alternative specification of a| is necessary. 
Here we investigate damping, alternative definitions of a| are left to a future study. 



Zhao and Mohseni [241 ] found in an a priori study using a dynamic proce- 
dure that a followed a linear variation from zero at the wall to a constant value of 
0.02 around ~ 10. In a later study [2^ they tested this distribution for a (now 
turning off the dynamic procedure and fixing a to follow the specified profile), 
but their results showed a significant overprediction of the spanwise velocity 
fiuctuations and also some overprediction of the vertical velocity fiuctuations 
(together suggesting high streamwise vorticity). Otherwise their results were not 
significantly better than a standard LES. This is not surprising because their 
a distribution was determined using an a priori study, which does not account 
for the feedback of the model on the flow. For example, in the previous section 
we saw that the effect of having the streamwise vortices closer to the wall is to 
increase skin friction. This is the type of effect you will not see in an a priori study. 

Following Zhao and Mohseni we used a linear variation of a and from 
zero at the wall to a constant value at a specifled wall-normal distance. Here the 
damping was applied over the region < 60 instead of the region y"*" < 10 that 
was used in their study. This choice was motivated by the study by Jimenez that 
demonstrated in the region y^ < 60 the low-speed streaks are a critical part of 
the autonomous cycle of near- wall turbulence [52]. Thus we consider this to be the 
'streak-affected' region, and since the problem is related to the velocity gradients 
from the streaks, du/dz, it is logical to apply the damping factor through this 
region. Other values were tested along with exponential damping profiles instead 
of the simple linear one. The shape of the damping profile was found to be 
insignificant, with the wall-normal distance being the important factor. Damping 
over the streak affected region was found to consistently provide the best results. 
The damping function used was 

= if .^<60 

•'^^ ^ \ 1 otherwise. ^ ' 

Results with damping are shown in Figure [8] for the isotropic and anisotropic 
models. These are now for full channels with (Lx,Ly,Lz) = (47r, 2, 27r/3) 
and {Nx,Ny,Nz) = (32,48,32). For the isotropic model was specified as 
oi^ = f (y^) (0.02)^, where / (y+) is given by equation ([65]) and 0.02 is the value of 
a away from the wall determined by Zhao and Mohseni [24] . We also tried a value 
of 0.04, which they determined for the i?eT-=180 channel from scaling arguments, 
but we found this was too high to yield reasonable results. 
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For the anisotropic model a\ = f (y+) Ch\ was used for and while 
was set to zero. Initially a C value of 1/6 was used such that it is a damped 
version of the case where the filter width is twice the grid spacing. However, 
with this value the logarithmic law was still underpredicted (in terms of the 
y- intercept), so the results reported here used C = 1/12. Even with this value 
the skin friction is still overpredicted. To have the same physical equivalent a for 
the isotropic and anisotropic models you would need to use C = 1/24 or A = /i. 
This was tested (results not shown) and it did bring the mean velocity profile 
into good agreement with the DNS data, but it does not seem to make good 
physical sense because this would mean the filter width is equal to the grid spacing. 

We can see in Figure [8] that damping removes the problem with the high 
spanwise fluctuations, and improves some quantities slightly (eg. mean flow pro- 
file, shear stress and streamwise velocity fluctuations) but overall the differences 
between the no model and NS-a model results are very small when damping is 
used. In Figure [8] we also show results from a simulation using the anisotropic 
Leray model (with djUj = enforced), that is the Aij + Bij terms in equation 
(I59p . The Leray model does not have the same vortex tilting properties that the 
NS-a model has. We can see by comparing the Leray model results in Figure [8] 
with those from the NS-a model in Figures [T] and [6] that the Leray model does 
not suffer the same underprediction of the mean flow that the NS-a model does, 
reinforcing the fact that it is the Cij term that is causing the problem. The tilting 
term, UkdiUk, which combines with the modifled pressure gradient in equation 
(??) to form the dj{Cij) term in the model, is the unique feature of the NS-a 
model. 

4.4.4. Helicity PDFs 

In the previous two sections we have seen that the NS-a model has a tendency to 
tilt vorticity close to the wall, and that with damping this impact is reduced. We 
now address the question of how the model changes the vorticity and velocity fields 



away from the wall. To do this we look at the relative helicity, defined as [53[] 



u ■ u , , 

/i = — n— r. 66 
\u\\uj\ 

The helicity is related to the nonlinear term u x u through the identity 

{u-uf {u X 

I. .191. .ic^ + = 1- (67) 



Thus by looking at the helicity, we can also examine the non-linearity of the model. 
Given that the NS-a equations are described as having a reduced nonlinearity [l3] , 
we expect they may also have high helicity. For the NS-a model in the LES-template 
we can write 

-+uxu = -v[, + -u.uy .V^u - H-^ j (68) 



and investigate the smoothed relative helicity 



June 13, 2010 17;26 Journal of Turbulence hpchanS 



20 

The focus now is on the region away from the wall, thus the PDFs shown in Figure 
[9] were measured at ~ 90. We can see for the minimal channel the NS-a model 
has two shoulders near ±1 indicating a higher probability of increased helicity 
(reduced nonlinearity) relative to the Navier-Stokes equations. Values of the mean- 
squared helicity, h"^, are given in Table [2j A uniform distribution would correspond 
to h"^ = 0.33. For the full channel we can see when no model is used the PDF is too 
peaked and h"^ is too low, but when the damped NS-a model is used the results 
are closer to the minimal channel DNS (which is itself close to the full channel 
DNS of Rogers and Moin [s^]). Because we expect the minimal channel DNS to 
be representative of a full channel at a finer mesh spacing, this suggests that the 
NS-a model can produce helicity statistics on a coarse mesh that are comparable 
to those from a finer mesh without a model. 

5. Conclusions 

In this paper the NS-a model has been investigated for a fully turbulent channel 
flow. To begin we derived the model using Hamilton's principle. Using this 
derivation it is straightforward to see how the model can be extended to different 
physical situations. For example, compressible flow or geophysical flows, how 
the model could be altered by using a higher order expansion in the definition 
of the velocity fluctuation, or different definitions of ^. The definition of a| 
used in practice should be consistent with that used in the derivation. For 
example, we found that in our application of the NS-a model to the channel 
flow, when a| is based on the mesh the a^ values are too large and this leads to 
excessive tilting of spanwise vorticity into the streamwise direction in the near 
wall region. This is because we derived the model assuming that a| follows an 
advection equation, and technically we should have solved an advection equation 
to determine a|, instead of basing it on the computational mesh. The fact 
that the magnitudes of the values used for a^ with damping are close the those 
one would obtain from solving an advection equation [s^] reinforces this statement. 

Given the significant impact the NS-a model has on the vorticity field, we 
feel that future studies should investigate how the NS-a model affects the resolved 
flow vortices (see da Silva et al. [HI for a study of this nature for other subgrid 
models). Such a study would have practical implications as well, for example in 
applications where it is the size, strength and location of the fluid vortices that is 
of interest. 
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Table 1. Mesh parameters for the minimal channel flow. The first column are 
the {x,y,z) dimensions non-dimensionalized by the channel half-height, the 
second column is the number of mesh points in each direction, and the last 
three are the channel dimensions in wall units (non-dimensionalized by the 
shear velocity, u-r, and the viscosity, ly). 
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hy (min) / hy (max) 
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2,0.37r) 


(16,64,16) 


35.3 


0.875/11.6 


10.6 
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2,0.37r) 


(24,64,24) 


23.6 


0.875/11.6 


7.07 


III 
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2,0.37r) 


(32,64,32) 


17.7 


0.875/11.6 


5.30 


IV 




2,0.37r) 


(32,128,32) 


17.7 


0.424/5.80 


5.30 
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Table 2. Values of the mean-squared helicity 

based on the smoothed velocity and vorticity, 
h^, for the DNS and for the NS-a model. 







model 




DNS (minimal channel, mcsli IV) 


0.31 


NS-a (minimal channel, mesh IV) 


0.37 


NS-a (damping, full channel) 


0.29 


no model (full channel) 


0.21 



;26 Journal of Turbulence hpchan5 



23 



20 



+ 




(a) Mean flow 



Figure 1. Mean flow profiles for different meshes with different values of the smoothing length scale a, in 
terms of parameter C from equation II60I I. Dash-dotted line, C = 1/6, (16,64, 16); dashed line, C = 1/6 , 
(24, 64, 24); dotted Une, C = 2/3, (32, 64, 32); soUd C = 1/6, (32, 128, 32); symbols full channel DNS [SOl . 
To look at the effect of refining the mesh while keeping the subgrid resolution constant the dash-dotted, 
dashed and solid lines can bo compared. To look at the effect of keeping the physical size of constant 
while refining the mesh (hence increasing the subgrid resolution) the dash-dotted and dotted lines can be 
compared. 
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Figure 2. Contours of the instantaneous spanwise velocity for the minimal channel for the (32, 128, 32) 
mesh. There is more small scale activity with the NS-a model. 
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Figure 3. Energy transfer TgQg = niijdjUi. Dashed line, T^q^^; solid line, T^Q^g] dash-dotted line 
'^SGSC' dotted line, total transfer, Tsgsa + Tsgsb ~ Tsgsc- See equation I I59I I for the definitions of the 
A, B and C terms. The peak of the total energy transfer here is correct in magnitude, but too close to the 
wall. 
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(a) Streamwise 




(b) Spanwise 



Figure 4. (a) Streamwise and (b) Spanwise vorticity profiles. Solid line, no model (24, 64, 24); dashed line, 
NS-o model (32,128,32); dotted line, NS-a model (24,64,24); symbols full channel DNS [g. The peak 
corresponding to the streamwise vortices in (a) is very close to the wall, j/+ 9, for the NS-« model. The 
spanwise vorticity is also overpredicted. 
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(a) DNS (b) NS-a 



Figure 5. Two point correlations of the streamwise velocity in the spanwise direction. Left panel is DNS 
and right panel is the NS-a model. Solid line is at j/+ « 7, dashed line is 2/+ f» 18. The streak spa<;ing, 
indicated by the minimum of the two point correlation, is much smaller for the NS-a model than for the 
DNS. 
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Figure 6. Streamwise velocity fluctuations. Dash-dotted line, C = 1/6, (16, 64, 16); dashed line, C = 1/6, 
(24, 64, 24); dotted line, C = 2/3, (32, 64, 32); solid line C = 1/6, (32, 128, 32); symbols full channel DNS 
[SOU . To look at the effect of refining the mesh while keeping the subgrid resolution the dash-dotted, 
dashed and solid lines can be compared. To look at the effect of keeping the physical size of o| constant 
while refining the mesh (hence increasing the subgrid resolution) the dash-dotted and dotted lines can be 
compared. 
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Figure 7. Strcamwisc vorticity PDFs. Left column is DNS and right column is NS-a model. Notice that 
the peak corresponding to streamwise vortices (at 20° is very prominent for the NS-ct model close to 
the wall at j/+ = 5.4, while it doesn't appear in the DNS until 2/+ = 18. 
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(e) Shear stress 



Figure 8. Mean flow, rms and shear stress profiles; no model (solid); anisotropic NS-a with damping 
(dashed); isotropic NS-a with damping (dash-dot); anisotropic Leray without damping (dotted). Symbols 
are DNS data [SOll . Note in particular that the mean flow proflle and skin friction for the damped isotropic 
model are in very good agreement with the DNS. 
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(a) minimal channel, 2/+ f» 90 (b) full channel, j/+ Ri 90 



Figure 9. PDFs of the relative helicity, (a) minimal channel; DNS (mesh IV), solid; NS-a (mesh IV), 
dashed; (b) full channel; NS-a with damping (anisotropic model), dashed; no model, dash-dotted;, solid 
line is from the minimal channel DNS, shown for comparison. In (b) we can see that the NS-a helicity 
PDF is very similar to that from the DNS, whereas when no model is used the shape of the helicity PDF 
is incorrect. 
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